#include "gtest/gtest.h"
#include "DataStructure/Common/Constant.h"
#include "../../IntegrationFormula/Gauss2DTria3P.h"
#include "../Tria6NodeShpFunc.h"

using SG::DataStructure::Common::TOLERANCE;
using SG::Algebra::Matrixd;
using SG::FEMSolver::Element::Tria6NodeShpFunc;
using SG::DataStructure::FEM::IntegPointArray;

// -------------------------------------------------------------
//    Kronecker-Delta Test
//    验证：在每个节点自身坐标处，对应形函数=1，其余=0
// -------------------------------------------------------------
TEST(Tria6NodeShpFuncTest, shpFuncKroneckerDeltaTest)
{
    // 六节点二次三角形标准单元局部坐标
    IntegPointArray nodes(6);
    nodes[0].m_x = 0.0;  nodes[0].m_y = 0.0;  // Node 1
    nodes[1].m_x = 1.0;  nodes[1].m_y = 0.0;  // Node 2
    nodes[2].m_x = 0.0;  nodes[2].m_y = 1.0;  // Node 3
    nodes[3].m_x = 0.5;  nodes[3].m_y = 0.0;  // Node 4 (1–2 midpoint)
    nodes[4].m_x = 0.5;  nodes[4].m_y = 0.5;  // Node 5 (2–3 midpoint)
    nodes[5].m_x = 0.0;  nodes[5].m_y = 0.5;  // Node 6 (3–1 midpoint)

    Tria6NodeShpFunc shpFunc;
    shpFunc.Compute(nodes);

    // 每个节点坐标处，对应形函数值=1，其余=0
    for (std::size_t i = 0; i < 6; ++i)
    {
        const auto& Ni = shpFunc.GetShapeFunction (static_cast<int>(i));
        for (std::size_t j = 0; j < 6; ++j)
        {
            double expected = (i == j) ? 1.0 : 0.0;
            ASSERT_NEAR(Ni(0, j), expected, TOLERANCE);
        }
    }
}

// -------------------------------------------------------------
//    Partition of Unity Test
//    验证：任意点处所有形函数之和 = 1
// -------------------------------------------------------------
TEST(Tria6NodeShpFuncTest, shpFuncPartitionOfUnityTest)
{
    // 采用三点高斯积分点测试,实际上取任一点均可
    SG::FEMSolver::Element::Gauss2DTria3P integCalc;
    auto gaussPts = integCalc.GetPoints();

    Tria6NodeShpFunc shpFunc;
    shpFunc.Compute(gaussPts);

    for (std::size_t i = 0; i < gaussPts.size(); ++i)
    {
        const auto& Ni = shpFunc.GetShapeFunction(static_cast<int>(i));
        double sum = 0.0;
        for (int j = 0; j < 6; ++j)
            sum += Ni(0, j);

        ASSERT_NEAR(sum, 1.0, TOLERANCE);
    }
}

// -------------------------------------------------------------
//    Derivative Consistency Test
//    验证：dN/dξ + dN/dη + dN/d(1−ξ−η) = 0（几何约束）
// -------------------------------------------------------------
TEST(Tria6NodeShpFuncTest, shpFuncDerivativeConsistencyTest)
{
    SG::FEMSolver::Element::Gauss2DTria3P integCalc;
    auto gaussPts = integCalc.GetPoints();

    Tria6NodeShpFunc shpFunc;
    shpFunc.Compute(gaussPts);

    for (std::size_t i = 0; i < gaussPts.size(); ++i)
    {
        const auto& dN = shpFunc.GetShapeFunctionDerivative(static_cast<int>(i));

        // 对6个形函数的导数，检查行和应为0（因为 L1+L2+L3=1 恒等式）
        double sumDNDxi = 0.0;
        double sumDNDdeta = 0.0;
        for (int j = 0; j < 6; ++j)
        {
            sumDNDxi  += dN(0, j);
            sumDNDdeta += dN(1, j);
        }

        ASSERT_NEAR(sumDNDxi, 0.0, TOLERANCE);
        ASSERT_NEAR(sumDNDdeta, 0.0, TOLERANCE);
    }
}